import data from 86 files .nc4¶
- date: 2019-01-01
- variables:
- latitudes
- longitudes
- CO2_std
- CO2_sdev
- dimensions:
- time
- lat
- lon
- attributes:
- year
- month
- day
- time
- variables:
Normailize data ???¶
# import libs
import os
import netCDF4 as nc4
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
import itertools
# Ruta al directorio que contiene los archivos netCDF
directorio_archivos = 'descargas'
# Lista de archivos netCDF en el directorio
archivos_netCDF = [os.path.join(directorio_archivos, nombre_archivo) for nombre_archivo in os.listdir(directorio_archivos) if nombre_archivo.endswith('.nc4')]
# create a dict de dimensiones [91,144]
row = np.arange(0,91)
col = np.arange(0,144)
dict_dim = dict()
def media_armonica(valores):
inversos = [1 / x for x in valores if x != 0] # Calcula los inversos de los valores no iguales a 0
if len(inversos) == 0:
return 0 # Si no hay valores válidos, devuelve 0
return len(inversos) / sum(inversos)
for i,j in itertools.product(row,col):
dict_dim[(i,j)] = {"lat": [], "lon": [], "co2_std": [], "co2_sdev": []}
# dict_dim[(i,j,"lon")] = []
# dict_dim[(i,j,"co2_std")] = []
# dict_dim[(i,j,"co2_sdev")] = []
for i in range(len(archivos_netCDF)):
# Abre el archivo netCDF
data = nc4.Dataset(archivos_netCDF[i],"r")
latitude = data.variables['Latitude'][:]
longitude = data.variables['Longitude'][:]
co2_std = data.variables['mole_fraction_of_carbon_dioxide_in_free_troposphere'][:]
co2_sdev = data.variables['mole_fraction_of_carbon_dioxide_in_free_troposphere_sdev'][:]
# print(latitude)
# print(longitude)
# print(co2_std)
# print(co2_sdev)
lon, lat = np.meshgrid(longitude, latitude)
# Aplana las mallas de coordenadas y la matriz de co2_free_troposphere a vectores unidimensionales
lon_flat = lon #.flatten()
lat_flat = lat #.flatten()
# create data.fream
df_lat = pd.DataFrame(lat_flat).fillna(0)
df_lon = pd.DataFrame(lon_flat).fillna(0)
df_co2_std = pd.DataFrame(co2_std).fillna(0)
df_co2_sdev = pd.DataFrame(co2_sdev).fillna(0)
for i,j in itertools.product(row,col):
aux_lat = df_lat.iloc[i,j]
dict_dim[(i,j)]["lat"].append(aux_lat)
dict_dim[(i,j)]["lon"].append(df_lon.iloc[i,j])
dict_dim[(i,j)]["co2_std"].append(df_co2_std.iloc[i,j])
dict_dim[(i,j)]["co2_sdev"].append(df_co2_sdev.iloc[i,j])
# Crear un data frame de los puntos unicos lat, lon para la dimension i,j
# print(dict_dim[(0,0)]["lat"])
columns = ["co2_std", "co2_sdev", "lat", "lon"]
df_dim = pd.DataFrame(columns=columns)
glat, glon = [], []
for i,j in itertools.product(row,col):
# print(dict_dim[(i,j)]["lat"])
# print(dict_dim[(i,j)]["lon"])
# print(dict_dim[(i,j)]["co2_std"])
# print(dict_dim[(i,j)]["co2_sdev"])
# print("==================================")
# print("==================================")
# print("==================================")
# print("=============================
x = media_armonica(dict_dim[(i,j)]["co2_sdev"])
y = media_armonica(dict_dim[(i,j)]["co2_std"])
#print(pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0])
glat.append(pd.Series(dict_dim[(i,j)]["lat"]).unique()[0]), glon.append(pd.Series(dict_dim[(i,j)]["lon"]).unique()[0])
min_lat_lon = min([len(pd.Series(dict_dim[(i,j)]["lat"]).unique()),len(pd.Series(dict_dim[(i,j)]["lon"]).unique())])
if min_lat_lon == 1:
df_dim.loc[len(glat)-1] = [x, y, pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]
# ({"co2_std": x, "co2_sdev": y, "lat": [pd.Series(dict_dim[(i,j)]["lat"]).unique()[0]], "lon": [pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]})
else:
df_dim.loc[len(glat)-1] = [x, y, pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]
print("no es unico".upper().ljust(30,"="))
# print()
# value.append(df_lat.iloc[i,j])
# value.append(df_lon.iloc[i,j])
# normalize df_dim
x = range(6)
np.std(x)
np.mean(x)
2.5
df_dim["co2_sdev"]
0 0.000395
1 0.000407
2 0.000401
3 0.000394
4 0.000403
...
13099 0.000000
13100 0.000000
13101 0.000000
13102 0.000000
13103 0.000000
Name: co2_sdev, Length: 13104, dtype: float64
#create a geodataframe with geopandas and shapely: lat, lon -> glat, glon
#from matplotlib.patches import Polygon
#geometry = [gdp.Polygon(xy) for xy in zip(glon, glat)]
geometry = [Point(xy) for xy in zip(glon, glat)]
crs = 'EPSG:4326'
gdf = gpd.GeoDataFrame(df_dim, crs=crs, geometry=geometry)
min(df_dim["co2_std"])
float(max(df_dim["co2_std"]))
df_dim.query("co2_std >= 0.000005")
| co2_std | co2_sdev | lat | lon | |
|---|---|---|---|---|
| 269 | 0.000005 | 0.000396 | 88.0 | 132.5 |
#describe the data
gdf.describe()
| co2_std | co2_sdev | lat | lon | |
|---|---|---|---|---|
| count | 13104.000000 | 13104.000000 | 13104.000000 | 13104.000000 |
| mean | 0.000002 | 0.000330 | -0.005495 | -1.250000 |
| std | 0.000001 | 0.000148 | 52.528319 | 103.924508 |
| min | 0.000000 | 0.000000 | -90.000000 | -180.000000 |
| 25% | 0.000002 | 0.000394 | -46.000000 | -90.625000 |
| 50% | 0.000003 | 0.000395 | 0.000000 | -1.250000 |
| 75% | 0.000003 | 0.000396 | 46.000000 | 88.125000 |
| max | 0.000005 | 0.000419 | 89.500000 | 177.500000 |
gdf.info()
<class 'geopandas.geodataframe.GeoDataFrame'> Index: 13104 entries, 0 to 13103 Data columns (total 5 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 co2_std 13104 non-null float64 1 co2_sdev 13104 non-null float64 2 lat 13104 non-null float64 3 lon 13104 non-null float64 4 geometry 13104 non-null geometry dtypes: float64(4), geometry(1) memory usage: 614.2 KB
import plotly.express as px
fig = px.histogram(df_dim, x="co2_std", nbins=80, title="Histograma de co2_std")
fig.show()
import plotly.express as px
fig = px.histogram(df_dim, x="co2_sdev", nbins=200, title="Histograma de co2_sdev")
fig.show()
Basado en el segso observado en los graficos, optamos por eliminar valores donde co2_std -> 0¶
gdf2 = gdf.query("co2_std >= 0.000000001") #.plot(column='co2_std', cmap='OrRd', figsize=(15, 15), legend=True, scheme='quantiles')
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.colors as mcolors
import matplotlib.cm as cm
gdf2["co2_std"]
8 0.000003
91 0.000003
119 0.000004
144 0.000002
145 0.000002
...
10939 0.000002
10940 0.000003
10941 0.000002
10942 0.000003
10943 0.000003
Name: co2_std, Length: 10803, dtype: float64
Distribucion de datos despues de eliminar valores iguales o cercanos a 0¶
import plotly.express as px
fig = px.histogram(gdf2, x="co2_std", nbins=80, title="Histograma de co2_std")
fig.show()
import plotly.express as px
fig = px.histogram(gdf2, x="co2_sdev", nbins=200, title="Histograma de co2_sdev")
fig.show()
Analisis de estacionariedad¶
lati = gdf2["lat"]
loni = gdf2["lon"]
latlon = list(zip(lati, loni))
latlon = [str(x) for x in latlon]
fig,ax = plt.subplots()
plt.plot(latlon,gdf2["co2_std"],color='red')
plt.xticks(range(0,len(latlon),(round(len(latlon)/20))))
plt.xticks(rotation=60)
plt.xlabel("(Latitud, Longitud)")
plt.ylabel("Valor Co2_std")
plt.title("CO2 troposferico")
plt.show()
lati = gdf2["lat"]
loni = gdf2["lon"]
latlon = list(zip(lati, loni))
latlon = [str(x) for x in latlon]
fig,ax = plt.subplots()
plt.plot(latlon,gdf2["co2_sdev"],color='purple')
plt.xticks(range(0,len(latlon),(round(len(latlon)/20))))
plt.xticks(rotation=60)
plt.xlabel("(Latitud, Longitud)")
plt.ylabel("Valor Co2_sdev")
plt.title("desviacion estandar CO2 troposferico")
plt.show()
Pruebas de estacionariedad¶
Para esto utilizaremos pruebas de raiz unitaria. Una raíz unitaria es una característica de los procesos que evolucionan a través del tiempo y que puede causar problemas en inferencia estadística en modelos de series de tiempo.
Un proceso estocástico lineal tiene una raíz unitaria si el valor de la raíz de la ecuación característica del proceso es igual a 1, por lo tanto tal proceso es no estacionario. Si las demás raíces de la ecuación característica se encuentran dentro del círculo unitario es decir, tienen un valor absoluto menor a uno. entonces la primera diferencia del proceso es estacionaria.
Algunas pruebas de raiz unitarias son:¶
Argumented dickey fuller¶
La prueba de Dickey-Fuller aumentada es una prueba de hipótesis. La hipótesis nula es que la serie de tiempo no es estacionaria, y la alternativa es que la serie es estacionaria. Por lo tanto, necesitamos encontrar un valor p lo suficientemente bajo como para rechazar nuestra hipótesis nula, lo que sugiere que la serie es estacionaria.
Phillips Perron¶
La prueba Phillips perron es una prueba de raiz unitaria basada en la prueba de dickey fuller y su misma hipotesis nula. y si el valor p es menor a la significanica rechazamos la hipotesis nula.
Kwiatkowski Phillips Schmidt Shin¶
Contrariamente a la mayoría de las pruebas de raíz unitaria , la presencia de una raíz unitaria no es la hipótesis nula sino la alternativa. Además, en la prueba KPSS, la ausencia de una raíz unitaria no es una prueba de estacionariedad sino, por diseño, de estacionariedad debil. Esta es una distinción importante ya que es posible que una serie de tiempo no sea estacionaria, no tenga raíz unitaria pero sea debilmente estacionaria. por lo tanto si el valor p es menor a la significancia se dice que la serie es debilmente estacionaria.
Dickey Fuller¶
La Prueba de Dickey-Fuller busca determinar la existencia o no de raíces unitarias en una serie de tiempo. La hipótesis nula de esta prueba es que existe una raíz unitaria en la serie. El Modelos de regresión puede ser escrito como:
$$ \nabla y_{t}=(\rho-1)y_{t-1}+u_{t}=\delta y_{t-1}+u_{t} $$
Donde ∇ es el operador de primera diferencia. Este modelo puede ser estimado y las pruebas para una raíz unitaria son equivalentes a pruebas δ = 0 (donde δ = ρ - 1). y al igual que en la prueba aumentada para rechazar la hipotesis nula el valor p obtenido de la prueba debe ser menor a nuestra significancia.
# Librerias
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from arch.unitroot import PhillipsPerron,KPSS,DFGLS,ADF
from prettytable import PrettyTable
# ignore future warnings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
gdf2.columns
Index(['co2_std', 'co2_sdev', 'lat', 'lon', 'geometry'], dtype='object')
Pruebas para CO2 Troposférico std¶
serie = gdf2.iloc[:,0]
adf = ADF(serie)
pp = PhillipsPerron(serie)
kpss = KPSS(serie)
df = DFGLS(serie)
x = PrettyTable()
x.field_names = ["Unitary root test","Test Statistic", "P-value", "Critical Value for 5%", "Alternative Hypothesis:"]
x.add_row(["Argumented Dickey Fuller", adf.stat, adf.pvalue, adf.critical_values["5%"], adf.alternative_hypothesis])
x.add_row(["Phillips Perron", pp.stat, pp.pvalue, pp.critical_values["5%"], pp.alternative_hypothesis])
x.add_row(["Kwiatkowski-Phillips-Schmidt-Shin", kpss.stat, kpss.pvalue, kpss.critical_values["5%"], kpss.alternative_hypothesis])
x.add_row(["Dickey-Fuller GLS", df.stat, df.pvalue, df.critical_values["5%"], df.alternative_hypothesis])
x
| Unitary root test | Test Statistic | P-value | Critical Value for 5% | Alternative Hypothesis: |
|---|---|---|---|---|
| Argumented Dickey Fuller | -9.27058995597038 | 1.3301353820184625e-15 | -2.8618085021011126 | The process is weakly stationary. |
| Phillips Perron | -101.07287313284054 | 0.0 | -2.861807607138271 | The process is weakly stationary. |
| Kwiatkowski-Phillips-Schmidt-Shin | 3.5542876734971465 | 0.0001 | 0.4614 | The process contains a unit root. |
| Dickey-Fuller GLS | -8.486833870023165 | 6.730213808921115e-14 | -1.9456491396867541 | The process is weakly stationary. |
En resumen, las pruebas de Dickey-Fuller Aumentada, Phillips-Perron y Dickey-Fuller GLS indican que el proceso es débilmente estacionario, ya que los valores P son muy cercanos a cero, y los estadísticos de prueba están por debajo de los valores críticos.
Por otro lado, la prueba KPSS indica que el proceso contiene una unidad raíz, lo que significa que no es estacionario. Esto crea una discrepancia entre las pruebas, y en este caso, se daría preferencia a las pruebas de Dickey-Fuller Aumentada, Phillips-Perron y Dickey-Fuller GLS, que indican estacionariedad. Sin embargo, es importante considerar la interpretación en el contexto de los datos y los objetivos del análisis. Esta incertidumbre debido a la discrepancia en los resultados, nos motivan a considerar otras técnicas de análisis.
Pruebas para la desviacion estandar del Co2¶
serie = gdf2.iloc[:,1]
adf = ADF(serie)
pp = PhillipsPerron(serie)
kpss = KPSS(serie)
df = DFGLS(serie)
x = PrettyTable()
x.field_names = ["Unitary root test","Test Statistic", "P-value", "Critical Value for 5%", "Alternative Hypothesis:"]
x.add_row(["Argumented Dickey Fuller", adf.stat, adf.pvalue, adf.critical_values["5%"], adf.alternative_hypothesis])
x.add_row(["Phillips Perron", pp.stat, pp.pvalue, pp.critical_values["5%"], pp.alternative_hypothesis])
x.add_row(["Kwiatkowski-Phillips-Schmidt-Shin", kpss.stat, kpss.pvalue, kpss.critical_values["5%"], kpss.alternative_hypothesis])
x.add_row(["Dickey-Fuller GLS", df.stat, df.pvalue, df.critical_values["5%"], df.alternative_hypothesis])
x
| Unitary root test | Test Statistic | P-value | Critical Value for 5% | Alternative Hypothesis: |
|---|---|---|---|---|
| Argumented Dickey Fuller | -4.438444574142901 | 0.00025362816868106694 | -2.8618085769516295 | The process is weakly stationary. |
| Phillips Perron | -25.853410981967464 | 0.0 | -2.861807607138271 | The process is weakly stationary. |
| Kwiatkowski-Phillips-Schmidt-Shin | 12.415548071750424 | 0.0001 | 0.4614 | The process contains a unit root. |
| Dickey-Fuller GLS | -0.7686066800043696 | 0.3939327633407306 | -1.945649700954751 | The process is weakly stationary. |
En este caso, la Prueba de Dickey-Fuller Aumentada y la Prueba de Phillips-Perron indican que el proceso es débilmente estacionario, ya que los valores P son muy bajos, y los estadísticos de prueba están por debajo de los valores críticos. Estas pruebas respaldan la estacionariedad de los datos.
Por otro lado, la Prueba de Kwiatkowski-Phillips-Schmidt-Shin (KPSS) sugiere que el proceso contiene una unidad raíz, lo que significa que no es estacionario. Esta prueba está en contradicción con las otras dos pruebas.
La Prueba de Dickey-Fuller GLS no es tan concluyente, ya que el valor P es relativamente alto, y el estadístico de prueba no rechaza la hipótesis de no estacionariedad, pero tampoco la confirma con firmeza.
En contexto¶
Dado el contexto, es importante considerar que la estacionariedad de estas series temporales puede ser crucial en análisis geoestadísticos relacionados con el CO2 atmosférico, ya que las tendencias y patrones temporales pueden influir en la interpretación de los datos y en la toma de decisiones.
Las pruebas de raíz unitaria que has realizado indican que la Prueba de Dickey-Fuller Aumentada y la Prueba de Phillips-Perron sugieren que las series de datos "co2_std" y "co2_sdev" son débilmente estacionarias, lo que implica que pueden ser utilizadas en análisis estadísticos con la suposición de estacionariedad.
Por otro lado, la Prueba de Kwiatkowski-Phillips-Schmidt-Shin (KPSS) sugiere que las series de datos contienen una unidad raíz, lo que implicaría no estacionariedad. Esta discrepancia entre las pruebas podría deberse a diferentes supuestos y enfoques de las pruebas.
En esta situación, consideramos las siguientes opciones:
Evaluar visualmente los datos: Realiza gráficos de las series temporales para identificar tendencias, estacionalidades u otros patrones visuales que puedan ayudarte a tomar una decisión más informada sobre la estacionariedad.
Realizar más pruebas: Podrías considerar realizar pruebas adicionales o explorar técnicas estadísticas más avanzadas para evaluar la estacionariedad de las series de datos.
Diferenciación: Si no puedes llegar a una conclusión clara sobre la estacionariedad de las series de datos, podrías aplicar técnicas de diferenciación para transformar los datos y hacerlos estacionarios antes de continuar con tu análisis geoestadístico.
Logaritmos: En nuestro caso particular y debido que los valores son pequeños, aplicar logaritmos aún puede ser una estrategia útil para estabilizar la varianza y reducir la magnitud de las fluctuaciones en la serie temporal. Sin embargo, es importante tener en cuenta que, en este caso, el logaritmo natural (ln) puede ser más apropiado, ya que tiende a funcionar mejor con valores pequeños. Como datos son muy pequeños, es probable que tengan una magnitud similar a la de los valores logaritmicos resultantes. Esto puede facilitar la interpretación de los resultados y hacer que los datos transformados sean más adecuados para análisis groestadísticos, incluyendo pruebas de estacionariedad.
kpss.alternative_hypothesis
'The process contains a unit root.'
kpss.summary()
| Test Statistic | 3.554 |
| P-value | 0.000 |
| Lags | 58 |
Trend: Constant
Critical Values: 0.74 (1%), 0.46 (5%), 0.35 (10%)
Null Hypothesis: The process is weakly stationary.
Alternative Hypothesis: The process contains a unit root.
# Crea una figura de Matplotlib
fig, ax = plt.subplots(figsize=(12, 8))
# Dibuja el mapa del mundo de GeoPandas como fondo
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.boundary.plot(ax=ax, linewidth=1, color='black')
# Utiliza Seaborn para dibujar los puntos y personalizar la paleta de colores
sns.scatterplot(ax=ax,x= 'lon', y= 'lat', data = gdf2, hue = 'co2_std', palette='coolwarm')
# Agrega opacidad a los puntos para diferenciar la concentración
alpha = 0.3
ax.collections[0].set_alpha(alpha)
# Personaliza el aspecto del mapa
ax.set_title('Distribución de CO2 en la Troposfera Libre con Opacidad')
plt.axis('off') # Desactiva los ejes
sm = cm.ScalarMappable(cmap='coolwarm')
sm.set_norm(mcolors.Normalize(vmin=min(gdf2["co2_std"]), vmax=max(gdf2["co2_std"])))
# Agrega una barra de colores (colorbar)
cbar = plt.colorbar(sm, ax=ax)
cbar.set_label('Concentración de CO2 (co2_std)')
# Muestra el mapa con los datos
plt.show()
/tmp/ipykernel_2109/3664182266.py:5: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
# Crea una figura de Matplotlib
fig, ax = plt.subplots(figsize=(12, 8))
# Dibuja el mapa del mundo de GeoPandas como fondo
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.boundary.plot(ax=ax, linewidth=1, color='black')
# Utiliza Seaborn para dibujar los puntos y personalizar la paleta de colores
sns.scatterplot(ax=ax,x= 'lon', y= 'lat', data = gdf2, hue = 'co2_sdev', palette='viridis')
# Agrega opacidad a los puntos para diferenciar la concentración
alpha = 0.3
ax.collections[0].set_alpha(alpha)
# Personaliza el aspecto del mapa
ax.set_title('Distribución de CO2_sdev en la Troposfera Libre con Opacidad')
plt.axis('off') # Desactiva los ejes
sm = cm.ScalarMappable(cmap='viridis')
sm.set_norm(mcolors.Normalize(vmin=min(gdf2["co2_sdev"]), vmax=max(gdf2["co2_sdev"])))
# Agrega una barra de colores (colorbar)
cbar = plt.colorbar(sm, ax=ax)
cbar.set_label('Concentración de CO2 (co2_sdev)')
# Muestra el mapa con los datos
plt.show()
/tmp/ipykernel_2109/1589977332.py:5: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
gdf2.columns
Index(['co2_std', 'co2_sdev', 'lat', 'lon', 'geometry'], dtype='object')
import numpy as np
import pandas as pd
from geostatsmodels import utilities, variograms, model, kriging, geoplot
# Extrae las indices y los valores de "co2_std" del GeoDataFrame
indices = gdf2.geometry.apply(lambda geom: [geom.x, geom.y]).to_list()
values = gdf2[["lon","lat" , "co2_std"]]
P = values.to_numpy()
P
array([[-1.60000000e+02, 8.95000000e+01, 3.03339266e-06],
[ 4.75000000e+01, 8.95000000e+01, 3.42648696e-06],
[ 1.17500000e+02, 8.95000000e+01, 3.83136512e-06],
...,
[ 1.72500000e+02, -6.00000000e+01, 2.16231664e-06],
[ 1.75000000e+02, -6.00000000e+01, 2.50185760e-06],
[ 1.77500000e+02, -6.00000000e+01, 2.77377723e-06]])
pw = utilities.pairwise(P)
geoplot.hscattergram(P,pw,10,5)
geoplot.hscattergram(P,pw,20,5)
geoplot.hscattergram(P,pw,30,5)
gdf["co2_std"].mean()
# Calcula el semivariograma
tolerance = gdf["co2_std"].std()
lags = np.arange( tolerance, gdf["co2_std"].mean(), tolerance * 2 )
sill = np.var(P[:,2])
svm = model.semivariance( model.spherical, [ 4000, sill ] )
geoplot.semivariogram( P, lags, tolerance, model=svm )
You have more than 10,000 data points, this might take a minute.
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) /workspaces/Geo/note.ipynb Celda 25 line 7 <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X30sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a> sill = np.var(P[:,2]) <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X30sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a> svm = model.semivariance( model.spherical, [ 4000, sill ] ) ----> <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X30sdnNjb2RlLXJlbW90ZQ%3D%3D?line=6'>7</a> geoplot.semivariogram( P, lags, tolerance, model=svm ) File ~/.python/current/lib/python3.10/site-packages/geostatsmodels/geoplot.py:87, in semivariogram(data, lags, tol, model) 75 ''' 76 Input: (data) NumPy array with three columns, the first two 77 columns should be the x and y coordinates, and (...) 84 Output: empirical semivariogram 85 ''' 86 # h, sv = variograms.semivariogram(data, lags, tol) ---> 87 vdata = variograms.semivariogram(data, lags, tol) 88 h, sv = vdata[0], vdata[1] 89 sill = np.var(data[:, 2]) File ~/.python/current/lib/python3.10/site-packages/geostatsmodels/variograms.py:71, in semivariogram(data, lags, tol) 63 def semivariogram(data, lags, tol): 64 ''' 65 Input: (data) NumPy array where the first two columns 66 are the spatial coordinates, x and y (...) 69 Output: (sv) <2xN> NumPy array of lags and semivariogram values 70 ''' ---> 71 return variogram(data, lags, tol, 'semivariogram') File ~/.python/current/lib/python3.10/site-packages/geostatsmodels/variograms.py:121, in variogram(data, lags, tol, method) 119 v = [covariance(data, indices) for indices in index] 120 # bundle the semivariogram values with their lags --> 121 return np.c_[lags, v].T File ~/.local/lib/python3.10/site-packages/numpy/lib/index_tricks.py:418, in AxisConcatenator.__getitem__(self, key) 414 # concatenate could do cast, but that can be overriden: 415 objs = [array(obj, copy=False, subok=True, 416 ndmin=ndmin, dtype=final_dtype) for obj in objs] --> 418 res = self.concatenate(tuple(objs), axis=axis) 420 if matrix: 421 oldndim = res.ndim ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 1 and the array at index 1 has size 0
# Ajusta un modelo de semivariograma
# Puedes elegir diferentes modelos (spherical, exponential, etc.) y ajustar sus parámetros.
vgm = model.vgm(gdf2, "co2_std", "exponential", 10, 10)
# Grafica el semivariograma y el modelo ajustado
import matplotlib.pyplot as plt
plt.scatter(svm['lag'], svm['val'], label='Semivariogram', s=60)
plt.plot(vgm['distance'], vgm['gamma'], label='Model', color='red')
plt.xlabel('Lag Distance')
plt.ylabel('Semivariance')
plt.legend()
plt.show()
--------------------------------------------------------------------------- NameError Traceback (most recent call last) /workspaces/Geo/note.ipynb Celda 26 line 3 <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a> # Ajusta un modelo de semivariograma <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=1'>2</a> # Puedes elegir diferentes modelos (spherical, exponential, etc.) y ajustar sus parámetros. ----> <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=2'>3</a> vgm = model.vgm(gdf2, "co2_std", "exponential", 10, 10) <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a> # Grafica el semivariograma y el modelo ajustado <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a> import matplotlib.pyplot as plt NameError: name 'model' is not defined